# A Rainbow Ceiling? Sexual Orientation and Party Leader Evaluations

#Cozza, Joseph Francesco, Gonzalo Di Landro, Andrea Aldrich, and Zeynep Somer-Topcu


# British Journal of political Science, 2025


library(haven)
library(tidyr)
library(plyr)
library(dplyr)
library(stringr)
library(ggplot2)
library(cregg)
library(knitr)
library(xtable)
library(gridExtra)
library(projoint)
library(stargazer)
library(vtable)
library(labelled)
library(ltm)
library(tools)
library(scales)  
library(patchwork)
library(epiDisplay)

setwd(dirname(rstudioapi::getActiveDocumentContext()$path))
rm(list = ls(all = TRUE))
load("final_data.RData")


# Figure 1 (full results) ---------------------------------------------------------------


# AMCE Estimation
h1_amce <- cj(
  final_data,
  subsleg_index ~ selectorate + candidates + polling + lsex + lorientation + lage + lexperience + lperformance,
  estimate = "amce",
  id = ~id
)

# Linear Model (for reference/stats)
mod <- lm(subsleg_index ~ selectorate + candidates + polling + lsex + lorientation + lage + lexperience + lperformance, data = final_data)
nrow(model.frame(mod))
table(final_data$subsleg_index, useNA = "always")

# AMCE Table
df_h1_amce <- data.frame(h1_amce)[, 3:ncol(h1_amce)]
names(df_h1_amce) <- toTitleCase(tolower(names(df_h1_amce)))
amce_table <- xtable(df_h1_amce)
print(amce_table, type = "latex", include.rownames = FALSE, caption.placement = "top")


# AMCE Plot (Filtered by 'lorientation')

amce_lorientation <- subset(h1_amce, feature == "lorientation")

p1_h1 <- ggplot(amce_lorientation, aes(x = estimate, y = level)) +
  geom_point(position = position_dodge(0.5), size = 2) +
  geom_errorbar(aes(xmin = lower, xmax = upper), position = position_dodge(0.5), width = 0.1, size = 1) +
  geom_vline(xintercept = 0, linetype = "dashed") +
  xlim(-0.5, 0.5) +
  theme_bw() +
  ggtitle("Average Marginal Component Effect") +
  labs(x = "Estimated AMCE", y = "") +
  theme(plot.title = element_text(size = 10))


# MM Estimation

h1_mm <- cj(
  final_data,
  subsleg_index ~ selectorate + candidates + polling + lsex + lorientation + lage + lexperience + lperformance,
  estimate = "mm",
  id = ~id
)

# MM Table

df_h1_mm <- data.frame(h1_mm)[, 3:ncol(h1_mm)]
names(df_h1_mm) <- toTitleCase(tolower(names(df_h1_mm)))
mm_table <- xtable(df_h1_mm)
print(mm_table, type = "latex", include.rownames = FALSE, caption.placement = "top")


# MM Plot (Filtered by 'lorientation')
mm_lorientation <- subset(h1_mm, feature == "lorientation")

p2_h1 <- ggplot(mm_lorientation, aes(x = estimate, y = level)) +
  geom_point(position = position_dodge(0.5), size = 2) +
  geom_errorbar(aes(xmin = lower, xmax = upper), position = position_dodge(0.5), width = 0.1, size = 1) +
  geom_vline(xintercept = 2.5, linetype = "dashed") +
  xlim(2, 3) +
  theme_bw() +
  ggtitle("Marginal Means") +
  labs(x = "Estimated MMs", y = "") +
  theme(plot.title = element_text(size = 10))

# Combined Plot
grid.arrange(p1_h1, p2_h1, ncol = 2)


# Figure 2 (full results) ----------------------------------------------------------------

# MM Estimation (By Sexual Orientation)
h1abc_mm <- cj(
  final_data,
  subsleg_index ~ lsex,
  id = ~id,
  estimate = "mm",
  by = ~lorientation
)

# Extract Upper Bounds for Plot Labels
gay_upper_woman      <- h1abc_mm$upper[h1abc_mm$level == "Woman" & h1abc_mm$lorientation == "Gay"]
straight_upper_woman <- h1abc_mm$upper[h1abc_mm$level == "Woman" & h1abc_mm$lorientation == "Straight"]
gay_upper_man        <- h1abc_mm$upper[h1abc_mm$level == "Man" & h1abc_mm$lorientation == "Gay"]
straight_upper_man   <- h1abc_mm$upper[h1abc_mm$level == "Man" & h1abc_mm$lorientation == "Straight"]

# MM Plot
plot_h1abc_mm <- ggplot(h1abc_mm, aes(x = estimate, y = level, color = lorientation)) +
  geom_point(position = position_dodge(0.5), size = 2) +
  geom_errorbar(aes(xmin = lower, xmax = upper),
                position = position_dodge(0.5), width = 0.2, size = 1) +
  geom_vline(xintercept = 2.5, linetype = "dashed") +
  annotate("text", x = gay_upper_woman + 0.02, y = 2.27, label = "Gay", color = "#00BFC4", hjust = 0.45, size = 3) +
  annotate("text", x = straight_upper_woman + 0.02, y = 2.05, label = "Straight", color = "#F8766D", hjust = 0.45, size = 3) +
  annotate("text", x = gay_upper_man + 0.02, y = 1.27, label = "Gay", color = "#00BFC4", hjust = 0.45, size = 3) +
  annotate("text", x = straight_upper_man + 0.02, y = 1.05, label = "Straight", color = "#F8766D", hjust = 0.45, size = 3) +
  scale_x_continuous(
    limits = c(2.00, 3.00),
    breaks = seq(2.00, 3.00, by = 0.25),
    labels = sprintf("%.2f", seq(2.00, 3.00, by = 0.25))
  ) +
  theme_bw() +
  ggtitle("Marginal Means") +
  labs(x = "Estimated MMs", y = "", color = "") +
  theme(plot.title = element_text(size = 10),
        legend.position = "none")

# Linear Model Summary (Optional)
mod <- lm(subsleg_index ~ lsex * lorientation, data = final_data)
nrow(model.frame(mod))
table(final_data$subsleg_index, useNA = "always")

# MM Table (LaTeX Output)
df_h1abc_mm <- data.frame(h1abc_mm)[, c(1, 4:(ncol(h1abc_mm) - 1))]
names(df_h1abc_mm) <- toTitleCase(tolower(names(df_h1abc_mm)))
mm_table <- xtable(df_h1abc_mm)
print(mm_table, type = "latex", include.rownames = FALSE, caption.placement = "top")

# MM Differences Estimation (Gay vs. Straight)
h1abc_diffmm <- cj(
  final_data,
  subsleg_index ~ lsex,
  id = ~id,
  estimate = "mm_differences",
  by = ~lorientation
)

# MM Differences Plot
plot_h1abc_diffmm <- ggplot(h1abc_diffmm, aes(x = estimate, y = level)) +
  geom_point(position = position_dodge(0.5), size = 2) +
  geom_errorbar(aes(xmin = lower, xmax = upper),
                position = position_dodge(0.5), width = 0.1, size = 1) +
  geom_vline(xintercept = 0, linetype = "dashed") +
  scale_x_continuous(
    limits = c(-0.5, 0.5),
    breaks = seq(-0.5, 0.5, by = 0.25),
    labels = sprintf("%.2f", seq(-0.5, 0.5, by = 0.25))
  ) +
  theme_bw() +
  ggtitle("Difference in Marginal Means") +
  labs(x = "MMs (Gay - Straight)", y = "") +
  theme(plot.title = element_text(size = 10),
        legend.position = "none")

# MM Differences Table (LaTeX Output)
df_h1abc_diffmm <- data.frame(h1abc_diffmm)[, c(1, 4:(ncol(h1abc_diffmm) - 1))]
names(df_h1abc_diffmm) <- toTitleCase(tolower(names(df_h1abc_diffmm)))
diffmm_table <- xtable(df_h1abc_diffmm)
print(diffmm_table, type = "latex", include.rownames = FALSE, caption.placement = "top")

# Combined Plot
grid.arrange(plot_h1abc_mm, plot_h1abc_diffmm, ncol = 2)

# Figure 3 (full Results) ----------------------------------------------------------------

# MMs Estimation
h2_mm <- cj(
  final_data,
  subsleg_index ~ lexperience,
  id = ~id,
  estimate = "mm",
  by = ~lorientation
)

# Extract upper bounds for each level by orientation
gay_upper_23      <- h2_mm$upper[h2_mm$level == "Member of Parliament for 23 years" & h2_mm$lorientation == "Gay"]
straight_upper_23 <- h2_mm$upper[h2_mm$level == "Member of Parliament for 23 years" & h2_mm$lorientation == "Straight"]

gay_upper_17      <- h2_mm$upper[h2_mm$level == "Member of Parliament for 17 years" & h2_mm$lorientation == "Gay"]
straight_upper_17 <- h2_mm$upper[h2_mm$level == "Member of Parliament for 17 years" & h2_mm$lorientation == "Straight"]

gay_upper_11      <- h2_mm$upper[h2_mm$level == "Member of Parliament for 11 years" & h2_mm$lorientation == "Gay"]
straight_upper_11 <- h2_mm$upper[h2_mm$level == "Member of Parliament for 11 years" & h2_mm$lorientation == "Straight"]

gay_upper_5       <- h2_mm$upper[h2_mm$level == "Member of Parliament for 5 years" & h2_mm$lorientation == "Gay"]
straight_upper_5  <- h2_mm$upper[h2_mm$level == "Member of Parliament for 5 years" & h2_mm$lorientation == "Straight"]

# Plot for MMs
plot_h2_mm <- ggplot(h2_mm, aes(x = estimate, y = level, color = lorientation)) +
  geom_point(position = position_dodge(0.5), size = 2) +
  geom_errorbar(aes(xmin = lower, xmax = upper), position = position_dodge(0.5), width = 0.2, size = 1) +
  
  # Labels for each experience level
  annotate("text", x = gay_upper_23 + 0.02, y = 4.3, label = "Gay", color = "#00BFC4", hjust = 0.45, size = 3) +
  annotate("text", x = straight_upper_23 + 0.02, y = 4.1, label = "Straight", color = "#F8766D", hjust = 0.45, size = 3) +
  
  annotate("text", x = gay_upper_17 + 0.02, y = 3.3, label = "Gay", color = "#00BFC4", hjust = 0.45, size = 3) +
  annotate("text", x = straight_upper_17 + 0.02, y = 3.1, label = "Straight", color = "#F8766D", hjust = 0.45, size = 3) +
  
  annotate("text", x = gay_upper_11 + 0.02, y = 2.3, label = "Gay", color = "#00BFC4", hjust = 0.95, size = 3) +
  annotate("text", x = straight_upper_11 + 0.02, y = 2.1, label = "Straight", color = "#F8766D", hjust = 0.45, size = 3) +
  
  annotate("text", x = gay_upper_5 + 0.02, y = 1.3, label = "Gay", color = "#00BFC4", hjust = 0.45, size = 3) +
  annotate("text", x = straight_upper_5 + 0.02, y = 1.1, label = "Straight", color = "#F8766D", hjust = 0.55, size = 3) +
  
  geom_vline(xintercept = 2.5, linetype = "dashed") +
  scale_x_continuous(
    limits = c(1.9, 3.1),
    breaks = seq(2.00, 3.00, by = 0.25),
    labels = sprintf("%.2f", seq(2.00, 3.00, by = 0.25))
  ) +
  theme_bw() +
  ggtitle("Marginal Means") +
  labs(y = "", x = "Estimated MMs", color = "") +
  theme(
    plot.title = element_text(size = 10),
    legend.position = "none"
  )

# Linear Model Summary
mod <- lm(subsleg_index ~ lexperience * lorientation, data = final_data)
nrow(model.frame(mod))
table(final_data$subsleg_index, useNA = "always")

# MMs Table (LaTeX)
df_h2_mm <- data.frame(h2_mm)[, c(1, 4:(ncol(h2_mm) - 1))]
names(df_h2_mm) <- toTitleCase(tolower(names(df_h2_mm)))
mm_table <- xtable(df_h2_mm)
print(mm_table, type = "latex", include.rownames = FALSE, caption.placement = "top")

# MM Differences Estimation
h2_diffmm <- cj(
  final_data,
  subsleg_index ~ lexperience,
  id = ~id,
  estimate = "mm_difference",
  by = ~lorientation
)

# Plot for Differences in MMs
plot_h2_diffmm <- ggplot(h2_diffmm, aes(x = estimate, y = level)) +
  geom_point(position = position_dodge(0.5), size = 2) +
  geom_errorbar(aes(xmin = lower, xmax = upper), position = position_dodge(0.5), width = 0.1, size = 1) +
  geom_vline(xintercept = 0, linetype = "dashed") +
  scale_x_continuous(
    limits = c(-0.50, 0.50),
    breaks = seq(-0.50, 0.50, by = 0.25),
    labels = sprintf("%.2f", seq(-0.50, 0.50, by = 0.25))
  ) +
  theme_bw() +
  ggtitle("Difference in Marginal Means") +
  labs(y = "", x = "MMs (Gay - Straight)") +
  theme(
    plot.title = element_text(size = 10),
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank()
  )

# MM Differences Table (LaTeX)
df_h2_diffmm <- data.frame(h2_diffmm)[, c(1, 4:(ncol(h2_diffmm) - 1))]
names(df_h2_diffmm) <- toTitleCase(tolower(names(df_h2_diffmm)))
diffmm_table <- xtable(df_h2_diffmm)
print(diffmm_table, type = "latex", include.rownames = FALSE, caption.placement = "top")

# Combined Plot
plot_h2_mm + plot_h2_diffmm + plot_layout(ncol = 2, widths = c(1, 1))


# Appendix B - Sample ---------------------------------------------------------------
rm(list = ls(all = TRUE))
load("final_data.RData")

descriptives <- unique(
  final_data[,c("id","gender","age","edlevel","edlevel_60_text",
                "lrself","region","most_important_char","bias_battery_1",
                "bias_battery_2","bias_battery_3","bias_battery_4",
                "bias_battery_5","bias_battery_6","group")])

#age
descriptives$age <- as.numeric(descriptives$age)
descriptives$age_cat <- NA
descriptives[which(descriptives$age>17&descriptives$age<30),"age_cat"] <- "18-29"
descriptives[which(descriptives$age>29&descriptives$age<50),"age_cat"] <- "30-49"
descriptives[which(descriptives$age>49),"age_cat"] <- ">49"
descriptives$age_cat <- factor(descriptives$age_cat, levels = c("18-29","30-49",">49"))

#bias battery
descriptives <- descriptives %>%
  mutate(
    across(
      all_of(c("bias_battery_1",
               "bias_battery_2",
               "bias_battery_3",
               "bias_battery_4",
               "bias_battery_5",
               "bias_battery_6")), 
      ~ifelse(. == "Disagree strongly", 1,
              ifelse(. == "Disagree", 2,
                     ifelse(. == "Neither agree nor disagree", 3,
                            ifelse(. == "Agree", 4,
                                   ifelse(. == "Agree strongly", 5, .)))))))


descriptives$bias_battery_1 <- as.numeric(descriptives$bias_battery_1)
descriptives$bias_battery_2 <- as.numeric(descriptives$bias_battery_2)
descriptives$bias_battery_3 <- as.numeric(descriptives$bias_battery_3)
descriptives$bias_battery_4 <- as.numeric(descriptives$bias_battery_4)
descriptives$bias_battery_5 <- as.numeric(descriptives$bias_battery_5)
descriptives$bias_battery_6 <- as.numeric(descriptives$bias_battery_6)

# now reverse coding those for which higher values mean more bias
descriptives <- descriptives %>%
  mutate(bias_battery_2 = 6 - bias_battery_2,
         bias_battery_5 = 6 - bias_battery_5)


#correcting age education (some respondents put the year they finished education instead of
# how old they were when they finished. We have to estimate how many years passed since they finished, and then subtract that number from their age)

#below the code to manually see visualize those 7 cases
tocorrect  <- descriptives[which(descriptives$edlevel_60_text>60),c("age","edlevel_60_text")]
tocorrect$edu <- tocorrect$age - (2023 - tocorrect$edlevel_60_text)

#correcting
descriptives$edlevel_60_text <- ifelse(descriptives$edlevel_60_text>60,
                                       descriptives$age - (2023 - descriptives$edlevel_60_text),
                                       descriptives$edlevel_60_text)

#coding 2 respondents as missing gender (sex)
descriptives[which(descriptives$gender=="In another way"),"gender"] <- NA
descriptives$woman <- ifelse(descriptives$gender=="Female",1,0)

descriptives$region <- factor(descriptives$region)

descriptives$young <- factor(ifelse(descriptives$age<50,"young","not-young"))




st(descriptives[,c("age","age_cat","woman","edlevel_60_text","lrself","region","bias_battery_1",
                   "bias_battery_2","bias_battery_3","bias_battery_4","bias_battery_5","bias_battery_6")], out="return")


#BES 2019 (last survey on a representative sample available)

bes <- haven::read_dta("bes_rps_2019_1.3.0.dta")
bes <- bes[,c("y09","Age","region")]

bes <- subset(bes,bes$y09>0&bes$y09<3)
bes$woman <- ifelse(bes$y09==2,1,0)

bes$age <- bes$Age
bes <- subset(bes,bes$age>0)

#age
bes$age_cat <- NA
bes[which(bes$age>17&bes$age<30),"age_cat"] <- "18-29"
bes[which(bes$age>29&bes$age<50),"age_cat"] <- "30-49"
bes[which(bes$age>49),"age_cat"] <- ">49"
bes$age_cat <- factor(bes$age_cat, levels = c("18-29","30-49",">49"))
bes$region <-  as_factor(bes$region)

st(bes[,c("age","age_cat","woman","region")],out = "return")










# Appendix C - Experimental manipulation -------------------------------------------
rm(list = ls(all = TRUE))
load("final_data.RData")

#Figure A1
plot((cj_freqs(
  final_data, ~ profile + selectorate+candidates + polling + lsex + lorientation + lage +        
    lexperience + lperformance, id=~id)))+ 
  ggplot2:: ggtitle("Distribution of Attribute Values")




# Appendix E.1 - Intra-Respondent Reliability estimation and correction ------------------------------------------------------
load("final_data.RData")
load("task5.RData")

# Getting outcome measures for task5
task5 <- task5[, c(1:3, 14:21)]

# Invert profile numbers to represent flipped profiles from task 1 for IRR estimation
task5$profile <- ifelse(task5$profile == 1, 2, 1)

# Save temporary version for merging with task 1
temp <- task5
temp$task <- 1

# Add repeated task outcomes
full <- merge(final_data, temp, by = c("id", "task", "profile"), all.x = TRUE)

# Check available repeated responses
nrow(full[which(!is.na(full$subsleg1_y)),]) # 9112
nrow(full[which(!is.na(full$subsleg2_y)),]) # 9088
nrow(full[which(!is.na(full$subsleg3_y)),]) # 9066
nrow(full[which(!is.na(full$subsleg4_y)),]) # 9116
nrow(full[which(!is.na(full$subsleg5_y)),]) # 9090

# IRR formula: IRR = 1 − 2τ(1 − τ)
# Example IRRs (from previous runs):
# subsleg 1: 0.683
# subsleg 2: 0.680
# subsleg 3: 0.670
# subsleg 4: 0.684
# subsleg 5: 0.697


#Subsleg item 1
projoint <- make_projoint_data(
  full,
  .attribute_vars = c("selectorate","candidates","polling","lsex","lorientation","lage","lexperience","lperformance"),
  .id_var = "id",
  .task_var = "task",
  .profile_var = "profile",
  .selected_var = "subsleg1_y",
  .selected_repeated_var = "rep_subsleg1_y",
  .fill = F
)
projoint <- read_labels(projoint, "own_shape_IRR_labs.csv")
output <- projoint(projoint, .remove_ties=F)
print(output)
summary(output)

#IRR design method
1- ((2 * output@tau) * (1-output@tau)) # design IRR = 0.6826836

# Plot: subsleg 1
plot(output, .estimates = "both") +
  theme_bw() +
  xlim(0.2, 0.8) +
  geom_vline(xintercept = 0.5) +
  labs(y = "", x = "Estimated MMs", color = "")

#Table: subsleg 1
tab <- data.frame(summary(output))
labs <- data.frame(output@labels)[, c(1, 2, 4)]
colnames(labs)[3] <- "att_level_choose"
tab <- merge(tab,labs, by = c("att_level_choose"), all = T)
tab <- tab[,2:ncol(tab)]
tab <- dplyr::select(tab, attribute, level, everything())
names(tab) <- toTitleCase(tolower(names(tab)))
tab$Estimand <- ifelse(tab$Estimand=="mm_uncorrected","MM uncorrected", "MM corrected")
colnames(tab)[5:7] <- c("Std.error", "Lower", "Upper")
tab <- xtable(tab) 
print(tab, type = "latex", include.rownames = FALSE, caption.placement = "top")





#Subsleg item 2
projoint <- make_projoint_data(
  full,
  .attribute_vars = c("selectorate","candidates","polling","lsex","lorientation","lage","lexperience","lperformance"),
  .id_var = "id",
  .task_var = "task",
  .profile_var = "profile",
  .selected_var = "subsleg2_y",
  .selected_repeated_var = "rep_subsleg2_y",
  .fill = F
)
projoint <- read_labels(projoint, "own_shape_IRR_labs.csv")
output <- projoint(projoint, .remove_ties=F)
print(output)
summary(output)


#IRR design method
1- ((2 * output@tau) * (1-output@tau)) # design IRR = 0.6802906

# Plot: subsleg 2
plot(output, .estimates = "both") +
  theme_bw() +
  xlim(0.2, 0.8) +
  geom_vline(xintercept = 0.5) +
  labs(y = "", x = "Estimated MMs", color = "")

#Table: subsleg 2
tab <- data.frame(summary(output))
labs <- data.frame(output@labels)[, c(1, 2, 4)]
colnames(labs)[3] <- "att_level_choose"
tab <- merge(tab,labs, by = c("att_level_choose"), all = T)
tab <- tab[,2:ncol(tab)]
tab <- dplyr::select(tab, attribute, level, everything())
names(tab) <- toTitleCase(tolower(names(tab)))
tab$Estimand <- ifelse(tab$Estimand=="mm_uncorrected","MM uncorrected", "MM corrected")
colnames(tab)[5:7] <- c("Std.error", "Lower", "Upper")
tab <- xtable(tab) 
print(tab, type = "latex", include.rownames = FALSE, caption.placement = "top")





#Subsleg item 3
projoint <- make_projoint_data(
  full,
  .attribute_vars = c("selectorate","candidates","polling","lsex","lorientation","lage","lexperience","lperformance"),
  .id_var = "id",
  .task_var = "task",
  .profile_var = "profile",
  .selected_var = "subsleg3_y",
  .selected_repeated_var = "rep_subsleg3_y",
  .fill = F
)
projoint <- read_labels(projoint, "own_shape_IRR_labs.csv")
output <- projoint(projoint, .remove_ties=F)
print(output)
summary(output)


#IRR design method
1- ((2 * output@tau) * (1-output@tau)) # design IRR = 0.6700091

# Plot: subsleg 3
plot(output, .estimates = "both") +
  theme_bw() +
  xlim(0.2, 0.8) +
  geom_vline(xintercept = 0.5) +
  labs(y = "", x = "Estimated MMs", color = "")


#Table: subsleg 3
tab <- data.frame(summary(output))
labs <- data.frame(output@labels)[, c(1, 2, 4)]
colnames(labs)[3] <- "att_level_choose"
tab <- merge(tab,labs, by = c("att_level_choose"), all = T)
tab <- tab[,2:ncol(tab)]
tab <- dplyr::select(tab, attribute, level, everything())
names(tab) <- toTitleCase(tolower(names(tab)))
tab$Estimand <- ifelse(tab$Estimand=="mm_uncorrected","MM uncorrected", "MM corrected")
colnames(tab)[5:7] <- c("Std.error", "Lower", "Upper")
tab <- xtable(tab) 
print(tab, type = "latex", include.rownames = FALSE, caption.placement = "top")





#Subsleg item 4
projoint <- make_projoint_data(
  full,
  .attribute_vars = c("selectorate","candidates","polling","lsex","lorientation","lage","lexperience","lperformance"),
  .id_var = "id",
  .task_var = "task",
  .profile_var = "profile",
  .selected_var = "subsleg4_y",
  .selected_repeated_var = "rep_subsleg4_y",
  .fill = F
)
projoint <- read_labels(projoint, "own_shape_IRR_labs.csv")
output <- projoint(projoint, .remove_ties=F)
print(output)
summary(output)

#IRR design method
1- ((2 * output@tau) * (1-output@tau)) # design IRR = 0.6838768

# Plot: subsleg 4
plot(output, .estimates = "both") +
  theme_bw() +
  xlim(0.2, 0.8) +
  geom_vline(xintercept = 0.5) +
  labs(y = "", x = "Estimated MMs", color = "")

#Table: subsleg 4
tab <- data.frame(summary(output))
labs <- data.frame(output@labels)[, c(1, 2, 4)]
colnames(labs)[3] <- "att_level_choose"
tab <- merge(tab,labs, by = c("att_level_choose"), all = T)
tab <- tab[,2:ncol(tab)]
tab <- dplyr::select(tab, attribute, level, everything())
names(tab) <- toTitleCase(tolower(names(tab)))
tab$Estimand <- ifelse(tab$Estimand=="mm_uncorrected","MM uncorrected", "MM corrected")
colnames(tab)[5:7] <- c("Std.error", "Lower", "Upper")
tab <- xtable(tab) 
print(tab, type = "latex", include.rownames = FALSE, caption.placement = "top")



#Subsleg item 5
projoint <- make_projoint_data(
  full,
  .attribute_vars = c("selectorate","candidates","polling","lsex","lorientation","lage","lexperience","lperformance"),
  .id_var = "id",
  .task_var = "task",
  .profile_var = "profile",
  .selected_var = "subsleg5_y",
  .selected_repeated_var = "rep_subsleg5_y",
  .fill = F
)
projoint <- read_labels(projoint, "own_shape_IRR_labs.csv")
output <- projoint(projoint, .remove_ties=F)
print(output)
summary(output)

#IRR design method
1- ((2 * output@tau) * (1-output@tau)) # design IRR = 0.6969973

# Plot: subsleg 5
plot(output, .estimates = "both") +
  theme_bw() +
  xlim(0.2, 0.8) +
  geom_vline(xintercept = 0.5) +
  labs(y = "", x = "Estimated MMs", color = "")

#Table: subsleg 5
tab <- data.frame(summary(output))
labs <- data.frame(output@labels)[, c(1, 2, 4)]
colnames(labs)[3] <- "att_level_choose"
tab <- merge(tab,labs, by = c("att_level_choose"), all = T)
tab <- tab[,2:ncol(tab)]
tab <- dplyr::select(tab, attribute, level, everything())
names(tab) <- toTitleCase(tolower(names(tab)))
tab$Estimand <- ifelse(tab$Estimand=="mm_uncorrected","MM uncorrected", "MM corrected")
colnames(tab)[5:7] <- c("Std.error", "Lower", "Upper")
tab <- xtable(tab) 
print(tab, type = "latex", include.rownames = FALSE, caption.placement = "top")








# Appendix E.2 - Validation of legitimacy index ----------------------------------------------------
rm(list = ls(all = TRUE))
load("final_data.RData")

# List of subsleg variables
subsleg_vars <- c("subsleg1_y", "subsleg2_y", "subsleg3_y", "subsleg4_y", "subsleg5_y")

# Corresponding labels
subsleg_labels <- c("Excluding Item 1:\nWhich of these leaders\nearned their position?",
                    "Excluding Item 2:\nWhich of these leaders\nwould be more effective\nin passing legislation?",
                    "Excluding Item 3:\nWhich of these leaders\nwill work harder\non behalf of their party?",
                    "Excluding Item 4:\nWhich of these leaders\nwould be more effective\nin unifying the party?",
                    "Excluding Item 5:\nWhich of these leaders\nwould help their party\nwin more seats\nin the next election?")

# Create a function to calculate the index excluding one variable
calc_index <- function(data, exclude_var = NULL) {
  if (is.null(exclude_var)) {
    # Use all variables
    included_vars <- subsleg_vars
  } else {
    # Exclude one variable
    included_vars <- subsleg_vars[subsleg_vars != exclude_var]
  }
  rowSums(data[included_vars])
}

# Store results
results <- data.frame()

# Calculate and store the model for the full index
final_data$subsleg_index <- calc_index(final_data)
full_model <- cj(
  final_data,
  subsleg_index ~ selectorate + candidates + polling + lsex + lorientation + lage + lexperience + lperformance,
  estimate = "amce",
  id = ~id
)
full_model_result <- full_model %>% 
  filter(feature == "lorientation") %>% 
  mutate(excluded_var = "Full Index")
results <- rbind(results, full_model_result)

# Loop through the combinations, excluding one variable at a time
for (i in seq_along(subsleg_vars)) {
  var <- subsleg_vars[i]
  label <- subsleg_labels[i]
  
  final_data$subsleg_index <- calc_index(final_data, var)
  
  model <- cj(
    final_data,
    subsleg_index ~ selectorate + candidates + polling + lsex + lorientation + lage + lexperience + lperformance,
    estimate = "amce",
    id = ~id
  )
  
  lorientation_result <- model %>% 
    filter(feature == "lorientation") %>% 
    mutate(excluded_var = label)
  
  results <- rbind(results, lorientation_result)
}

# Plot the results
ggplot(results, aes(x = factor(excluded_var, levels = c("Full Index", subsleg_labels)), y = estimate, color = ifelse(excluded_var == "Full Index", "Full Index", "Excluded Item"))) +
  geom_point(size = 2) +
  geom_errorbar(aes(ymin = lower, ymax = upper), width = 0.2, size = 1) + 
  labs(title = "",
       x = "",
       y = "Estimated AMCEs",
       color = "") +
  scale_color_manual(values = c("Full Index" = "black", "Excluded Item" = "#DF536B")) +
  theme_bw() +
  ylim(-0.5, 0.5) +
  geom_hline(yintercept = 0) +
  theme(legend.position = "none")+
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1))

# Cronbach's alpha 
cronbach_test <- subset(final_data, !is.na(final_data$subsleg1_y)&
                          !is.na(final_data$subsleg2_y)&
                          !is.na(final_data$subsleg3_y)&
                          !is.na(final_data$subsleg4_y)&
                          !is.na(final_data$subsleg5_y))
cronbach_test <- cronbach.alpha(
  cronbach_test[,c("subsleg1_y","subsleg2_y","subsleg3_y","subsleg4_y","subsleg5_y")])

cronbach_test
#Items: 5
#Sample units: 9050
#alpha: 0.844



# Mean and SD for each subsleg item and Index
library(knitr)

# List of variables to summarize
variables <- c("subsleg_index", "subsleg1_y", "subsleg2_y", "subsleg3_y", "subsleg4_y", "subsleg5_y")

# Create a function to compute summary statistics for a single variable
compute_summary <- function(data, variable) {
  data %>%
    summarize(
      Variable = variable,
      Mean = mean(.data[[variable]], na.rm = TRUE),
      Median = median(.data[[variable]], na.rm = TRUE),
      SD = sd(.data[[variable]], na.rm = TRUE),
      Min = min(.data[[variable]], na.rm = TRUE),
      Max = max(.data[[variable]], na.rm = TRUE),
      Missing = sum(is.na(.data[[variable]]))
    )
}

# Compute summaries for all variables
summary_table <- bind_rows(lapply(variables, function(var) compute_summary(final_data, var)))

# Display the summary table
summary_table

# Display the table in RStudio
kable(summary_table, format = "latex", caption = "Summary Statistics for Selected Variables", digits = 2, booktabs = TRUE)


# OLS Regression Models of Leader’s Sexual Orientation
library(sandwich)
library(lmtest)

# First Model: Base model with clustering
base_model <- lm(subsleg_index ~ lorientation, data = final_data)
clustered_se_base <- vcovCL(base_model, cluster = ~ id)
base_model_summary <- coeftest(base_model, vcov = clustered_se_base)
print(base_model_summary)

# Second Model: Fixed effects model
fe_model <- lm(subsleg_index ~ lorientation + as.factor(id), data = final_data)

# Adjust for clustering at the 'id' level in the FE model
clustered_se_fe <- vcovCL(fe_model, cluster = ~ id)

# Summarize the FE model with cluster-robust standard errors
fe_model_summary <- coeftest(fe_model, vcov = clustered_se_fe)
(fe_model_summary)[2,]


library(stargazer)
# Extract cluster-robust standard errors
base_model_se <- sqrt(diag(clustered_se_base))
fe_model_se <- sqrt(diag(clustered_se_fe))

# Create a Stargazer table with both models
stargazer(
  base_model, fe_model,
  se = list(base_model_se, fe_model_se), # Add robust standard errors
  title = "OLS Regression Models with Cluster-Robust Standard Errors",
  dep.var.labels = "Leader Evaluation Index",
  omit = "id",
  column.labels = c("Base Model", "Fixed Effects Model"),
  star.cutoffs = c(0.05, 0.01, 0.001),
  omit.stat=c("ser","f","adj.rsq"),
  add.lines = list(
    c("Respondent-clustered SE", "Yes","Yes"),
    c("Respondent FE", "No","Yes")),
  single.row = F,
  type="latex"
)


# Item missingness
missing_summary <- data.frame(
  Variable = c("subsleg1_y", "subsleg2_y", "subsleg3_y", "subsleg4_y", "subsleg5_y", "subsleg_index"),
  Non_Missing = c(
    sum(!is.na(final_data$subsleg1_y)),
    sum(!is.na(final_data$subsleg2_y)),
    sum(!is.na(final_data$subsleg3_y)),
    sum(!is.na(final_data$subsleg4_y)),
    sum(!is.na(final_data$subsleg5_y)),
    sum(!is.na(final_data$subsleg_index))
  ),
  Missing = c(
    sum(is.na(final_data$subsleg1_y)),
    sum(is.na(final_data$subsleg2_y)),
    sum(is.na(final_data$subsleg3_y)),
    sum(is.na(final_data$subsleg4_y)),
    sum(is.na(final_data$subsleg5_y)),
    sum(is.na(final_data$subsleg_index))
  )
)

missing_summary




# Appendix E.3 - Respondents who answered “The Position of Party Leader" ---------------------------------------------------
rm(list = ls(all = TRUE)) 
load("final_data.RData")

#missigness per value
final_data %>%
  group_by(mancheck1) %>%
  summarise(
    unique_ids = n_distinct(id),       
    total_obs = n()                   
  )

data_leaders_candidates <- final_data[!(final_data$mancheck1 %in% c("I don't know", "None of the above")), ]
data_leaders_candidates$passed_check <- ifelse(data_leaders_candidates$mancheck1=="The position of party leader",1,0)
data_leaders_candidates$passed_check <- factor(data_leaders_candidates$passed_check)

final_data$passed_check <- ifelse(final_data$mancheck1=="The position of party leader",1,0)
final_data$passed_check <- factor(final_data$passed_check)

#table
mod <- lm(subsleg_index ~ lorientation*passed_check, data=final_data)
nrow(model.frame(mod))

mod2 <- lm(subsleg_index ~ lorientation*passed_check, data=data_leaders_candidates)
nrow(model.frame(mod2))




# Results only for those who passed the check 
data_mancheck <- subset(final_data,final_data$mancheck1=="The position of party leader")
length(unique(data_mancheck$id)) # 605 respondents  


#Figure A10
# AMCE
h1_amce <- cj(
  data_mancheck,
  subsleg_index ~ selectorate + candidates + polling + lsex + lorientation + lage+ lexperience + lperformance,
  estimate = "amce",
  id = ~id)

#table
mod <- lm(subsleg_index ~ selectorate + candidates + polling + lsex + lorientation + lage+ lexperience + lperformance, data=data_mancheck)
nrow(model.frame(mod))
table(data_mancheck$subsleg_index,useNA = "always")
df_h1_amce <- data.frame(h1_amce)
df_h1_amce <- df_h1_amce[,3:ncol(df_h1_amce)]
names(df_h1_amce) <- toTitleCase(tolower(names(df_h1_amce)))
df_h1_amce <- xtable(df_h1_amce) #latex table
print(df_h1_amce, type = "latex", include.rownames = FALSE, 
      caption.placement = "top")


#plot
h1_amce <- h1_amce[which(h1_amce$feature=="lorientation"),]
p1_h1 <- ggplot(h1_amce,
                aes(x = estimate,
                    y = level)) +
  geom_point(position = position_dodge(width = 0.5),
             size = 2) +
  geom_errorbar(aes(xmin = lower, xmax = upper),
                position = position_dodge(width = 0.5),
                width = 0.2, size = 1) + 
  geom_vline(xintercept = 0) +
  xlim(-0.5,0.5)+
  theme_bw() +
  ggtitle("AMCEs") +
  xlab(NULL) +
  labs(
    y = "",   
    x="Estimated AMCEs",
    color = "")

# MM
h1_mm <- cj(
  data_mancheck,
  subsleg_index ~ selectorate + candidates + polling + lsex + lorientation + lage+ lexperience + lperformance,
  estimate = "mm",
  id = ~id)

#table
df_h1_mm <- data.frame(h1_mm)
df_h1_mm <- df_h1_mm[,3:ncol(df_h1_mm)]
names(df_h1_mm) <- toTitleCase(tolower(names(df_h1_mm)))
df_h1_mm <- xtable(df_h1_mm) #latex table
print(df_h1_mm, type = "latex", include.rownames = FALSE, 
      caption.placement = "top")

#plot
h1_mm <- h1_mm[which(h1_mm$feature=="lorientation"),]
p2_h1 <- ggplot(h1_mm,
                aes(x = estimate,
                    y = level)) +
  geom_point(position = position_dodge(width = 0.5),
             size = 2) +
  geom_errorbar(aes(xmin = lower, xmax = upper),
                position = position_dodge(width = 0.5),
                width = 0.2, size = 1) + 
  geom_vline(xintercept = 2.5) +
  xlim(2,3)+
  theme_bw() +
  ggtitle("MMs") +
  xlab(NULL) +
  labs(
    y = "",   
    x="Estimated MMs",
    color = "")


# combined plot (size = 3 x 6)
grid.arrange(p1_h1,p2_h1, ncol = 2)



#Figure A11
# MM
h1abc_mm <- cj(data_mancheck,
               subsleg_index ~ lsex,
               id = ~id, estimate = "mm",
               by = ~lorientation)

#table
mod <- lm(subsleg_index ~  lsex * lorientation, data=data_mancheck)
nrow(model.frame(mod))
table(data_mancheck$subsleg_index,useNA = "always")
df_h1abc_mm <- data.frame(h1abc_mm)
col <- ncol(df_h1abc_mm)-1
df_h1abc_mm <- df_h1abc_mm[,c(1,4:col)]
names(df_h1abc_mm) <- toTitleCase(tolower(names(df_h1abc_mm)))
df_h1abc_mm <- xtable(df_h1abc_mm) #latex table
print(df_h1abc_mm, type = "latex", include.rownames = FALSE, 
      caption.placement = "top")

#plot
plot_h1abc_mm <- ggplot(h1abc_mm,
                        aes(x = estimate,
                            y = level,
                            color = lorientation)) +
  geom_point(position = position_dodge(width = 0.5),
             size = 2) +
  geom_errorbar(aes(xmin = lower, xmax = upper),
                position = position_dodge(width = 0.5),
                width = 0.2, size = 1) + 
  geom_vline(xintercept = 2.5) +
  scale_x_continuous(breaks = c(2, 2.5, 3), limits = c(2, 3))+
  theme_bw() +
  ggtitle("MMs") +
  labs(y = "",
       x = "Estimated MMs",
       color = "") +
  theme(legend.position = "none") 

# Diff in MM
#table
h1abc_diffmm <- cj(data_mancheck,
                   subsleg_index ~ lsex,
                   id = ~id, estimate = "mm_differences",
                   by = ~lorientation)

df_h1abc_diffmm <- data.frame(h1abc_diffmm)
col <- ncol(df_h1abc_diffmm)-1
df_h1abc_diffmm <- df_h1abc_diffmm[,c(1,4:col)]
names(df_h1abc_diffmm) <- toTitleCase(tolower(names(df_h1abc_diffmm)))
df_h1abc_diffmm <- xtable(df_h1abc_diffmm) #latex table
print(df_h1abc_diffmm, type = "latex", include.rownames = FALSE, 
      caption.placement = "top")


#plot
plot_h1abc_diffmm <- ggplot(h1abc_diffmm,
                            aes(x = estimate,
                                y = level)) +
  geom_point(position = position_dodge(width = 0.5),
             size = 2) +
  geom_errorbar(aes(xmin = lower, xmax = upper),
                position = position_dodge(width = 0.5),
                width = 0.2, size = 1) + 
  geom_vline(xintercept = 0) +
  scale_x_continuous(breaks = c(-0.5, 0, 0.5), limits = c(-0.5,0.5))+
  theme_bw() +
  ggtitle("Difference in MMs") +
  xlab(NULL) +
  labs(
    y = "",
    x="MMs (Gay - Straight)")
# (size = 3 x 6)
grid.arrange(plot_h1abc_mm,plot_h1abc_diffmm, ncol=2)



#Figure A12

#MM 
#table
h2_mm <- cj(data_mancheck, subsleg_index~ lexperience, id = ~id, estimate = "mm", by = ~lorientation)
mod <- lm(subsleg_index ~ lexperience * lorientation, data=data_mancheck)
nrow(model.frame(mod))
table(final_data$subsleg_index,useNA = "always")
df_h2_mm <- data.frame(h2_mm)
col <- ncol(df_h2_mm)-1
df_h2_mm<- df_h2_mm[,c(1,4:col)]
names(df_h2_mm) <- toTitleCase(tolower(names(df_h2_mm)))
df_h2_mm <- xtable(df_h2_mm) #latex table
print(df_h2_mm, type = "latex", include.rownames = FALSE, 
      caption.placement = "top")

#plot
plot_h2_mm <- ggplot(h2_mm,
                     aes(x = estimate,
                         y = level,
                         color=lorientation)) +
  geom_point(position = position_dodge(width = 0.5),
             size = 2) +
  geom_errorbar(aes(xmin = lower, xmax = upper),
                position = position_dodge(width = 0.5),
                width = 0.2, size = 1) + 
  geom_vline(xintercept = 2.5) +
  scale_x_continuous(breaks = c(2, 2.5, 3), limits = c(1.75, 3.25))+
  theme_bw() +
  ggtitle("MMs") +
  xlab(NULL) +
  labs(
    y = "",   
    x="Estimated MMs",
    color = "")+
  theme(legend.position = "none") 


# MM diff
#table
h2_diffmm <- cj(data_mancheck, subsleg_index~ lexperience, id = ~id, estimate = "mm_difference", by = ~lorientation)
df_h2_diffmm <- data.frame(h2_diffmm)
col <- ncol(df_h2_diffmm)-1
df_h2_diffmm<- df_h2_diffmm[,c(1,4:col)]
names(df_h2_diffmm) <- toTitleCase(tolower(names(df_h2_diffmm)))
df_h2_diffmm <- xtable(df_h2_diffmm) #latex table
print(df_h2_diffmm, type = "latex", include.rownames = FALSE, 
      caption.placement = "top")

#plot
plot_h2_diffmm <- ggplot(h2_diffmm,
                         aes(x = estimate,
                             y = level)) +
  geom_point(position = position_dodge(width = 0.5),
             size = 2) +
  geom_errorbar(aes(xmin = lower, xmax = upper),
                position = position_dodge(width = 0.5),
                width = 0.2, size = 1) + 
  geom_vline(xintercept = 0) +
  scale_x_continuous(breaks = c(-0.5, 0, 0.5), limits = c(-0.5,0.5))+
  theme_bw() +
  ggtitle("Difference in MMs") +
  xlab(NULL) +
  labs(
    y = "",   
    x="MMs (Gay - Straight)",
    color = "")
# (size = 3 x 9)
grid.arrange(plot_h2_mm,plot_h2_diffmm, ncol=2)




#Figure A13
check_all_2 <- cj(
  final_data,
  subsleg_index ~ lorientation ,
  estimate = "mm_differences",
  by = ~passed_check,
  id = ~id
)

# Plot for `check_all`
p_check_all <- ggplot(check_all_2, aes(x = estimate, y = level)) +
  geom_point(position = position_dodge(width = 0.5), size = 2) +
  geom_errorbar(aes(xmin = lower, xmax = upper), width = 0.2, size = 1) + 
  geom_vline(xintercept = 0) +
  xlim(-0.5, 0.5) +
  theme_bw() +
  xlab(NULL) +
  labs(y = "", x = "Difference in MMs (passed check - did not pass)", color = "") 


















# Parties with queer leaders are perceived to be more left-leaning than parties with straight leaders

# AMCE
lr_amce <- cj(
  data_mancheck,
  lr ~ selectorate + candidates + polling + lsex + lorientation + lage+ lexperience + lperformance,
  id = ~id,
  estimate = "amce")

mod <- lm(lr ~ selectorate + candidates + polling + lsex + lorientation + lage+ lexperience + lperformance, data=data_mancheck)
nrow(model.frame(mod))
table(data_mancheck$lr,useNA = "always")

df_lr_amce <- data.frame(lr_amce)
df_lr_amce<- df_lr_amce[,3:ncol(df_lr_amce)]

names(df_lr_amce) <- toTitleCase(tolower(names(df_lr_amce)))
df_lr_amce <- xtable(df_lr_amce) #latex table

print(df_lr_amce, type = "latex", include.rownames = FALSE, 
      caption.placement = "top")




lr_amce <- lr_amce[which(lr_amce$feature=="lorientation"),]
p1_lr <- ggplot(lr_amce,
                aes(x = estimate,
                    y = level)) +
  geom_point(position = position_dodge(width = 0.5),
             size = 2) +
  geom_errorbar(aes(xmin = lower, xmax = upper),
                position = position_dodge(width = 0.5),
                width = 0.2, size = 1) + 
  geom_vline(xintercept = 0) +
  xlim(-0.8,0.8)+
  theme_bw() +
  theme(legend.position = "none") +
  ggtitle("AMCEs") +
  xlab(NULL) +
  labs(
    y = "",   
    x="Estimated AMCEs",
    color = "")

# MM
lr_mm <- cj(
  data_mancheck,
  lr ~ selectorate + candidates + polling + lsex + lorientation + lage+ lexperience + lperformance,
  id = ~id,
  estimate = "mm")

df_lr_mm <- data.frame(lr_mm)
df_lr_mm<- df_lr_mm[,3:ncol(df_lr_mm)]

names(df_lr_mm) <- toTitleCase(tolower(names(df_lr_mm)))
df_lr_mm <- xtable(df_lr_mm) #latex table

print(df_lr_mm, type = "latex", include.rownames = FALSE, 
      caption.placement = "top")


lr_mm <- lr_mm[which(lr_mm$feature=="lorientation"),]
p2_lr <- ggplot(lr_mm,
                aes(x = estimate,
                    y = level)) +
  geom_point(position = position_dodge(width = 0.5),
             size = 2) +
  geom_errorbar(aes(xmin = lower, xmax = upper),
                position = position_dodge(width = 0.5),
                width = 0.2, size = 1) + 
  geom_vline(xintercept = 5) +
  xlim(4,6)+
  theme_bw() +
  theme(legend.position = "none") +
  ggtitle("MMs") +
  xlab(NULL) +
  labs(
    y = "",   
    x="Estimated MMs",
    color = "")


# 3x6
grid.arrange(p1_lr,p2_lr, ncol = 2)
























# Appendix E.4 - Treatment effects conditional on respondents' bias towards gay politicians--------
rm(list = ls(all = TRUE))
load("final_data.RData")


descriptives <- unique(
  final_data[,c("id","gender","age","edlevel","edlevel_60_text",
                "lrself","region","most_important_char","bias_battery_1",
                "bias_battery_2","bias_battery_3","bias_battery_4",
                "bias_battery_5","bias_battery_6","group")])

#age
descriptives$age <- as.numeric(descriptives$age)
descriptives$age_cat <- NA
descriptives[which(descriptives$age>17&descriptives$age<30),"age_cat"] <- "18-29"
descriptives[which(descriptives$age>29&descriptives$age<50),"age_cat"] <- "30-49"
descriptives[which(descriptives$age>49),"age_cat"] <- ">49"
descriptives$age_cat <- factor(descriptives$age_cat, levels = c("18-29","30-49",">49"))

#bias battery
descriptives <- descriptives %>%
  mutate(
    across(
      all_of(c("bias_battery_1",
               "bias_battery_2",
               "bias_battery_3",
               "bias_battery_4",
               "bias_battery_5",
               "bias_battery_6")), 
      ~ifelse(. == "Disagree strongly", 1,
              ifelse(. == "Disagree", 2,
                     ifelse(. == "Neither agree nor disagree", 3,
                            ifelse(. == "Agree", 4,
                                   ifelse(. == "Agree strongly", 5, .)))))))


descriptives$bias_battery_1 <- as.numeric(descriptives$bias_battery_1)
descriptives$bias_battery_2 <- as.numeric(descriptives$bias_battery_2)
descriptives$bias_battery_3 <- as.numeric(descriptives$bias_battery_3)
descriptives$bias_battery_4 <- as.numeric(descriptives$bias_battery_4)
descriptives$bias_battery_5 <- as.numeric(descriptives$bias_battery_5)
descriptives$bias_battery_6 <- as.numeric(descriptives$bias_battery_6)


# now reverse coding those for which higher values mean more bias
descriptives <- descriptives %>%
  mutate(bias_battery_2 = 6 - bias_battery_2,
         bias_battery_5 = 6 - bias_battery_5)


#bias battery  
tab1(descriptives$group, sort.group = "decreasing", cum.percent = F, graph=F)


#Table A3
bias1 <- lm(bias_battery_1 ~group, data=descriptives)
bias2 <- lm(bias_battery_2 ~group, data=descriptives)
bias3 <- lm(bias_battery_3 ~group, data=descriptives)
bias4 <- lm(bias_battery_4 ~group, data=descriptives)
bias5 <- lm(bias_battery_5 ~group, data=descriptives)
bias6 <- lm(bias_battery_6 ~group, data=descriptives)

stargazer(bias1,bias2,bias3,bias4,bias5,bias6,
          title = "Answers to bias questions across bias-battery groups",
          dep.var.caption = "",
          dep.var.labels = c(
            "Hope female prime minister",
            "Men more capable",
            "Parties should more female MPs",
            "Hope gay prime minister",
            "Straight more capable",
            "Parties should more gay MPs"),
          star.cutoffs = c(0.05, 0.01, 0.001),
          omit.stat=c("ser","f"),
          type="text")



# coding Bias battery so that higher values mean more gender equal, and lower values more bias
final_data <- final_data %>%
  mutate(
    across(
      all_of(c("bias_battery_1",
               "bias_battery_2",
               "bias_battery_3",
               "bias_battery_4",
               "bias_battery_5",
               "bias_battery_6")), 
      ~ifelse(. == "Disagree strongly", 1,
              ifelse(. == "Disagree", 2,
                     ifelse(. == "Neither agree nor disagree", 3,
                            ifelse(. == "Agree", 4,
                                   ifelse(. == "Agree strongly", 5, .)))))))

final_data$bias_battery_1 <- as.numeric(final_data$bias_battery_1)
final_data$bias_battery_2 <- as.numeric(final_data$bias_battery_2)
final_data$bias_battery_3 <- as.numeric(final_data$bias_battery_3)
final_data$bias_battery_4 <- as.numeric(final_data$bias_battery_4)
final_data$bias_battery_5 <- as.numeric(final_data$bias_battery_5)
final_data$bias_battery_6 <- as.numeric(final_data$bias_battery_6)


# now reverse coding those for which higher values mean more bias
final_data <- final_data %>%
  mutate(bias_battery_2 = 6 - bias_battery_2,
         bias_battery_5 = 6 - bias_battery_5)

#creating an index for homophobia
final_data$homophobic_views_numeric <- final_data$bias_battery_4 + final_data$bias_battery_5 + final_data$bias_battery_6

table(final_data$homophobic_views_numeric)
median(final_data$homophobic_views_numeric,na.rm=T) # Median = 10

final_data$homophobic_views <- ifelse(final_data$homophobic_views_numeric < 10 ,
                                      "High bias against queer politicians",
                                      "Low bias against queer politicians")
final_data$homophobic_views <- as.factor(final_data$homophobic_views)


# MM
bias_mm <- cj(final_data,
              subsleg_index ~ lorientation,
              id = ~id, estimate = "mm",
              by = ~homophobic_views)


mod <- lm(subsleg_index ~ lorientation *homophobic_views, data=final_data)
nrow(model.frame(mod))
nrow(final_data[which(!is.na(final_data$subsleg_index)&!is.na(final_data$homophobic_views)),]) #8962
df_bias_mm <- data.frame(bias_mm)
col <- ncol(df_bias_mm)-1
df_bias_mm<- df_bias_mm[,c(1,4:col)]

names(df_bias_mm) <- toTitleCase(tolower(names(df_bias_mm)))
df_bias_mm <- xtable(df_bias_mm) #latex table
print(df_bias_mm, type = "latex", include.rownames = FALSE, 
      caption.placement = "top")


bias_mm$BY <- ifelse(bias_mm$BY == "Low bias against queer politicians", 
                     "Respondent with\nlow bias against\nqueer politicians", 
                     "Respondent with\nhigh bias against\nqueer politicians")


f1_homophobic <- ggplot(bias_mm,  aes(x = estimate, y = BY, color = level)) +
  geom_point(position = position_dodge(width = 0.5), size = 2) +
  geom_errorbar(aes(xmin = lower, xmax = upper),
                position = position_dodge(width = 0.5), width = 0.2, size = 1) + 
  geom_vline(xintercept = 2.5) +
  scale_x_continuous(breaks = c(2, 2.5, 3), limits = c(1.75, 3.25)) +
  theme_bw() +
  ggtitle("MMs") +
  theme(legend.position = "none", plot.title = element_text(size = 12)) + 
  labs(y = "", x = "Estimated MMs", color = "")



#differences in MMs
bias_diffmm <- cj(final_data,
                  subsleg_index ~ homophobic_views,
                  id = ~id, estimate = "mm_differences",
                  by = ~lorientation)


df_bias_diffmm <- data.frame(bias_diffmm)
col <- ncol(df_bias_diffmm)-1
df_bias_diffmm<- df_bias_diffmm[,c(1,4:col)]
names(df_bias_diffmm) <- toTitleCase(tolower(names(df_bias_diffmm)))
df_bias_diffmm <- xtable(df_bias_diffmm) #latex table
print(df_bias_diffmm, type = "latex", include.rownames = FALSE, 
      caption.placement = "top")



bias_diffmm$level <- ifelse(bias_diffmm$level == "Low bias against queer politicians", 
                            "Respondent with\nlow bias against\nqueer politicians", 
                            "Respondent with\nhigh bias against\nqueer politicians")


f2_homophobic <- ggplot(bias_diffmm,
                        aes(x = estimate,
                            y = level)) +
  geom_point(position = position_dodge(width = 0.5),
             size = 2) +
  geom_errorbar(aes(xmin = lower, xmax = upper),
                position = position_dodge(width = 0.5),
                width = 0.2, size = 1) + 
  geom_vline(xintercept = 0) +
  theme_bw() +
  scale_x_continuous(breaks = c(-1.5, 0, 1.5), limits = c(-1.5,1.5))+
  ggtitle("Difference in MMs") +
  theme(legend.position = "none",
        plot.title = element_text(size = 12)) + 
  xlab(NULL) +
  labs(
    y = "",
    x="MMs (Gay - Straight)")

grid.arrange(f1_homophobic,f2_homophobic,ncol=2)




# Respondents who passed the manipulation check
data_mancheck <- subset(final_data,final_data$mancheck1=="The position of party leader")


#coding Bias battery so that higher values mean more gender equal, and lower values more bias
data_mancheck <- data_mancheck %>%
  mutate(
    across(
      all_of(c("bias_battery_1",
               "bias_battery_2",
               "bias_battery_3",
               "bias_battery_4",
               "bias_battery_5",
               "bias_battery_6")), 
      ~ifelse(. == "Disagree strongly", 1,
              ifelse(. == "Disagree", 2,
                     ifelse(. == "Neither agree nor disagree", 3,
                            ifelse(. == "Agree", 4,
                                   ifelse(. == "Agree strongly", 5, .)))))))

data_mancheck$bias_battery_1 <- as.numeric(data_mancheck$bias_battery_1)
data_mancheck$bias_battery_2 <- as.numeric(data_mancheck$bias_battery_2)
data_mancheck$bias_battery_3 <- as.numeric(data_mancheck$bias_battery_3)
data_mancheck$bias_battery_4 <- as.numeric(data_mancheck$bias_battery_4)
data_mancheck$bias_battery_5 <- as.numeric(data_mancheck$bias_battery_5)
data_mancheck$bias_battery_6 <- as.numeric(data_mancheck$bias_battery_6)


# now reverse coding those for which higher values mean more bias
data_mancheck <- data_mancheck %>%
  mutate(bias_battery_2 = 6 - bias_battery_2,
         bias_battery_5 = 6 - bias_battery_5)

#creating an index for homophobia
data_mancheck$homophobic_views_numeric <- data_mancheck$bias_battery_4 + data_mancheck$bias_battery_5 + data_mancheck$bias_battery_6

table(data_mancheck$homophobic_views_numeric)
median(data_mancheck$homophobic_views_numeric,na.rm=T) # Median = 10

data_mancheck$homophobic_views <- ifelse(data_mancheck$homophobic_views_numeric < 10 ,
                                         "High bias against queer politicians",
                                         "Low bias against queer politicians")
data_mancheck$homophobic_views <- as.factor(data_mancheck$homophobic_views)


# MM
bias_mm <- cj(data_mancheck,
              subsleg_index ~ lorientation,
              id = ~id, estimate = "mm",
              by = ~homophobic_views)


mod <- lm(subsleg_index ~ lorientation *homophobic_views, data=data_mancheck)
nrow(model.frame(mod))
nrow(data_mancheck[which(!is.na(data_mancheck$subsleg_index)&!is.na(data_mancheck$homophobic_views)),]) #4748
df_bias_mm <- data.frame(bias_mm)
col <- ncol(df_bias_mm)-1
df_bias_mm<- df_bias_mm[,c(1,4:col)]

names(df_bias_mm) <- toTitleCase(tolower(names(df_bias_mm)))
df_bias_mm <- xtable(df_bias_mm) #latex table
print(df_bias_mm, type = "latex", include.rownames = FALSE, 
      caption.placement = "top")


bias_mm$BY <- ifelse(bias_mm$BY == "Low bias against queer politicians", 
                     "Respondent with\nlow bias against\nqueer politicians", 
                     "Respondent with\nhigh bias against\nqueer politicians")


f1_homophobic <- ggplot(bias_mm,  aes(x = estimate, y = BY, color = level)) +
  geom_point(position = position_dodge(width = 0.5), size = 2) +
  geom_errorbar(aes(xmin = lower, xmax = upper),
                position = position_dodge(width = 0.5), width = 0.2, size = 1) + 
  geom_vline(xintercept = 2.5) +
  scale_x_continuous(breaks = c(2, 2.5, 3), limits = c(1.75, 3.25)) +
  theme_bw() +
  ggtitle("MMs") +
  theme(legend.position = "none", plot.title = element_text(size = 12)) + 
  labs(y = "", x = "Estimated MMs", color = "")



#differences in MMs
bias_diffmm <- cj(data_mancheck,
                  subsleg_index ~ homophobic_views,
                  id = ~id, estimate = "mm_differences",
                  by = ~lorientation)


df_bias_diffmm <- data.frame(bias_diffmm)
col <- ncol(df_bias_diffmm)-1
df_bias_diffmm<- df_bias_diffmm[,c(1,4:col)]
names(df_bias_diffmm) <- toTitleCase(tolower(names(df_bias_diffmm)))
df_bias_diffmm <- xtable(df_bias_diffmm) #latex table
print(df_bias_diffmm, type = "latex", include.rownames = FALSE, 
      caption.placement = "top")



bias_diffmm$level <- ifelse(bias_diffmm$level == "Low bias against queer politicians", 
                            "Respondent with\nlow bias against\nqueer politicians", 
                            "Respondent with\nhigh bias against\nqueer politicians")


f2_homophobic <- ggplot(bias_diffmm,
                        aes(x = estimate,
                            y = level)) +
  geom_point(position = position_dodge(width = 0.5),
             size = 2) +
  geom_errorbar(aes(xmin = lower, xmax = upper),
                position = position_dodge(width = 0.5),
                width = 0.2, size = 1) + 
  geom_vline(xintercept = 0) +
  theme_bw() +
  scale_x_continuous(breaks = c(-1.5, 0, 1.5), limits = c(-1.5,1.5))+
  ggtitle("Difference in MMs") +
  theme(legend.position = "none",
        plot.title = element_text(size = 12)) + 
  xlab(NULL) +
  labs(
    y = "",
    x="MMs (Gay - Straight)")

#3x7
grid.arrange(f1_homophobic,f2_homophobic,ncol=2)





# Appendix E.5 - Treatment effects conditional on respondents' ideology---------------------------------------------------------------------
rm(list = ls(all = TRUE))
load("final_data.RData")

#creating left-wing indicator
median(final_data$lrself,na.rm = T) # Median = 5
final_data$right_participant <- factor(ifelse(final_data$lrself>5,"More conservative\nrespondent","Less conservative\nrespondent"))

#plot
rightwing_mm <- cj(final_data,
                   subsleg_index ~ lorientation,
                   id = ~id, estimate = "mm",
                   by = ~right_participant)


rightwing_mm$right_participant <- factor(rightwing_mm$right_participant, levels = c("More conservative\nrespondent","Less conservative\nrespondent"))

f1_rightwing <- ggplot(rightwing_mm,
                       aes(x = estimate,
                           y = right_participant,
                           color = level)) +
  geom_point(position = position_dodge(width = 0.5),
             size = 2) +
  geom_errorbar(aes(xmin = lower, xmax = upper),
                position = position_dodge(width = 0.5),
                width = 0.2, size = 1) + 
  geom_vline(xintercept = 2.5) +
  scale_x_continuous(breaks = c(2, 2.5, 3), limits = c(1.75, 3.25))+
  theme_bw() +
  ggtitle("MMs") +
  theme(legend.position = "none",
        plot.title = element_text(size = 12)) + 
  labs(
    y = "",   
    x="Estimated MMs",
    color = "")

#table
mod <- lm(subsleg_index ~ lorientation *right_participant, data=final_data)
nrow(model.frame(mod))
nrow(final_data[which(!is.na(final_data$subsleg_index)&!is.na(final_data$right_participant)),]) #9032
df_rightwing_mm <- data.frame(rightwing_mm)
col <- ncol(df_rightwing_mm)-1
df_rightwing_mm<- df_rightwing_mm[,c(1,4:col)]
names(df_rightwing_mm) <- toTitleCase(tolower(names(df_rightwing_mm)))
df_rightwing_mm <- xtable(df_rightwing_mm) #latex table
print(df_rightwing_mm, type = "latex", include.rownames = FALSE, 
      caption.placement = "top")

# Differences in MMs
rightwing_diffmm <- cj(final_data,
                       subsleg_index ~ right_participant,
                       id = ~id, estimate = "mm_differences",
                       by = ~lorientation)

#plot
rightwing_diffmm$level <- factor(rightwing_diffmm$level, levels=c("More conservative\nrespondent","Less conservative\nrespondent"))
f2_rightwing <- ggplot(rightwing_diffmm,
                       aes(x = estimate,
                           y = level)) +
  geom_point(position = position_dodge(width = 0.5), size = 2) +
  geom_errorbar(aes(xmin = lower, xmax = upper),
                position = position_dodge(width = 0.5),width = 0.2, size = 1) + 
  geom_vline(xintercept = 0) +
  theme_bw() +
  scale_x_continuous(breaks = c(-1.5, 0, 1.5), limits = c(-1.5,1.5))+
  ggtitle("Difference in MMs") +
  theme(legend.position = "none",
        plot.title = element_text(size = 12)) + 
  xlab(NULL) +
  labs(
    y = "",
    x="MMs (Gay - Straight)")

#table
df_rightwing_diffmm <- data.frame(rightwing_diffmm)
col <- ncol(df_rightwing_diffmm)-1
df_rightwing_diffmm<- df_rightwing_diffmm[,c(1,4:col)]
names(df_rightwing_diffmm) <- toTitleCase(tolower(names(df_rightwing_diffmm)))
df_rightwing_diffmm <- xtable(df_rightwing_diffmm) #latex table
print(df_rightwing_diffmm, type = "latex", include.rownames = FALSE, 
      caption.placement = "top")


# 3 x6
grid.arrange(f1_rightwing,f2_rightwing,ncol=2)



# Respondents who passed the manipulation check
data_mancheck <- subset(final_data,final_data$mancheck1=="The position of party leader")

#creating left-wing indicator
median(data_mancheck$lrself,na.rm = T) # Median = 5
data_mancheck$right_participant <- factor(ifelse(data_mancheck$lrself>5,"More conservative\nrespondent","Less conservative\nrespondent"))

#plot
rightwing_mm <- cj(data_mancheck,
                   subsleg_index ~ lorientation,
                   id = ~id, estimate = "mm",
                   by = ~right_participant)


rightwing_mm$right_participant <- factor(rightwing_mm$right_participant, levels = c("More conservative\nrespondent","Less conservative\nrespondent"))

f1_rightwing <- ggplot(rightwing_mm,
                       aes(x = estimate,
                           y = right_participant,
                           color = level)) +
  geom_point(position = position_dodge(width = 0.5),
             size = 2) +
  geom_errorbar(aes(xmin = lower, xmax = upper),
                position = position_dodge(width = 0.5),
                width = 0.2, size = 1) + 
  geom_vline(xintercept = 2.5) +
  scale_x_continuous(breaks = c(2, 2.5, 3), limits = c(1.75, 3.25))+
  theme_bw() +
  ggtitle("MMs") +
  theme(legend.position = "none",
        plot.title = element_text(size = 12)) + 
  labs(
    y = "",   
    x="Estimated MMs",
    color = "")

#table
mod <- lm(subsleg_index ~ lorientation *right_participant, data=data_mancheck)
nrow(model.frame(mod))
nrow(data_mancheck[which(!is.na(data_mancheck$subsleg_index)&!is.na(data_mancheck$right_participant)),]) #9032
df_rightwing_mm <- data.frame(rightwing_mm)
col <- ncol(df_rightwing_mm)-1
df_rightwing_mm<- df_rightwing_mm[,c(1,4:col)]
names(df_rightwing_mm) <- toTitleCase(tolower(names(df_rightwing_mm)))
df_rightwing_mm <- xtable(df_rightwing_mm) #latex table
print(df_rightwing_mm, type = "latex", include.rownames = FALSE, 
      caption.placement = "top")



# Differences in MMs
rightwing_diffmm <- cj(data_mancheck,
                       subsleg_index ~ right_participant,
                       id = ~id, estimate = "mm_differences",
                       by = ~lorientation)

#plot
rightwing_diffmm$level <- factor(rightwing_diffmm$level, levels=c("More conservative\nrespondent","Less conservative\nrespondent"))
f2_rightwing <- ggplot(rightwing_diffmm,
                       aes(x = estimate,
                           y = level)) +
  geom_point(position = position_dodge(width = 0.5), size = 2) +
  geom_errorbar(aes(xmin = lower, xmax = upper),
                position = position_dodge(width = 0.5),width = 0.2, size = 1) + 
  geom_vline(xintercept = 0) +
  theme_bw() +
  scale_x_continuous(breaks = c(-1.5, 0, 1.5), limits = c(-1.5,1.5))+
  ggtitle("Difference in MMs") +
  theme(legend.position = "none",
        plot.title = element_text(size = 12)) + 
  xlab(NULL) +
  labs(
    y = "",
    x="MMs (Gay - Straight)")

#table
df_rightwing_diffmm <- data.frame(rightwing_diffmm)
col <- ncol(df_rightwing_diffmm)-1
df_rightwing_diffmm<- df_rightwing_diffmm[,c(1,4:col)]
names(df_rightwing_diffmm) <- toTitleCase(tolower(names(df_rightwing_diffmm)))
df_rightwing_diffmm <- xtable(df_rightwing_diffmm) #latex table
print(df_rightwing_diffmm, type = "latex", include.rownames = FALSE, 
      caption.placement = "top")


# 3 x7
grid.arrange(f1_rightwing,f2_rightwing,ncol=2)





# Appendix E.6 - Treatment effects on ideological party positioning ---------------------------------

# AMCE
lr_amce <- cj(
  final_data,
  lr ~ selectorate + candidates + polling + lsex + lorientation + lage+ lexperience + lperformance,
  id = ~id,
  estimate = "amce")

mod <- lm(lr ~ selectorate + candidates + polling + lsex + lorientation + lage+ lexperience + lperformance, data=final_data)
nrow(model.frame(mod))
table(final_data$lr,useNA = "always")
df_lr_amce <- data.frame(lr_amce)
df_lr_amce<- df_lr_amce[,3:ncol(df_lr_amce)]
names(df_lr_amce) <- toTitleCase(tolower(names(df_lr_amce)))
df_lr_amce <- xtable(df_lr_amce) #latex table

print(df_lr_amce, type = "latex", include.rownames = FALSE, 
      caption.placement = "top")




lr_amce <- lr_amce[which(lr_amce$feature=="lorientation"),]
p1_lr <- ggplot(lr_amce,
                aes(x = estimate,
                    y = level)) +
  geom_point(position = position_dodge(width = 0.5),
             size = 2) +
  geom_errorbar(aes(xmin = lower, xmax = upper),
                position = position_dodge(width = 0.5),
                width = 0.2, size = 1) + 
  geom_vline(xintercept = 0) +
  xlim(-0.5,0.5)+
  theme_bw() +
  theme(legend.position = "none") +
  ggtitle("AMCEs") +
  xlab(NULL) +
  labs(
    y = "",   
    x="Estimated AMCEs",
    color = "")



# MM
lr_mm <- cj(
  final_data,
  lr ~ selectorate + candidates + polling + lsex + lorientation + lage+ lexperience + lperformance,
  id = ~id,
  estimate = "mm")


df_lr_mm <- data.frame(lr_mm)
df_lr_mm<- df_lr_mm[,3:ncol(df_lr_mm)]

names(df_lr_mm) <- toTitleCase(tolower(names(df_lr_mm)))
df_lr_mm <- xtable(df_lr_mm) #latex table

print(df_lr_mm, type = "latex", include.rownames = FALSE, 
      caption.placement = "top")


lr_mm <- lr_mm[which(lr_mm$feature=="lorientation"),]
p2_lr <- ggplot(lr_mm,
                aes(x = estimate,
                    y = level)) +
  geom_point(position = position_dodge(width = 0.5),
             size = 2) +
  geom_errorbar(aes(xmin = lower, xmax = upper),
                position = position_dodge(width = 0.5),
                width = 0.2, size = 1) + 
  geom_vline(xintercept = 5) +
  xlim(4,6)+
  theme_bw() +
  theme(legend.position = "none") +
  ggtitle("MMs") +
  xlab(NULL) +
  labs(
    y = "", 
    x="Estimated MMs",
    color = "")

# 3x6
grid.arrange(p1_lr,p2_lr, ncol = 2)



# Respondents who passed the manipulation check
data_mancheck <- subset(final_data,final_data$mancheck1=="The position of party leader")


lr_amce <- cj(
  data_mancheck,
  lr ~ selectorate + candidates + polling + lsex + lorientation + lage+ lexperience + lperformance,
  id = ~id,
  estimate = "amce")

mod <- lm(lr ~ selectorate + candidates + polling + lsex + lorientation + lage+ lexperience + lperformance, data=final_data)
nrow(model.frame(mod))
table(final_data$lr,useNA = "always")
df_lr_amce <- data.frame(lr_amce)
df_lr_amce<- df_lr_amce[,3:ncol(df_lr_amce)]
names(df_lr_amce) <- toTitleCase(tolower(names(df_lr_amce)))
df_lr_amce <- xtable(df_lr_amce) #latex table
print(df_lr_amce, type = "latex", include.rownames = FALSE, 
      caption.placement = "top")




lr_amce <- lr_amce[which(lr_amce$feature=="lorientation"),]
p1_lr <- ggplot(lr_amce,
                aes(x = estimate,
                    y = level)) +
  geom_point(position = position_dodge(width = 0.5),
             size = 2) +
  geom_errorbar(aes(xmin = lower, xmax = upper),
                position = position_dodge(width = 0.5),
                width = 0.2, size = 1) + 
  geom_vline(xintercept = 0) +
  xlim(-0.7,0.7)+
  theme_bw() +
  theme(legend.position = "none") +
  ggtitle("AMCEs") +
  xlab(NULL) +
  labs(
    y = "",   
    x="Estimated AMCEs",
    color = "")



# MM
lr_mm <- cj(
  data_mancheck,
  lr ~ selectorate + candidates + polling + lsex + lorientation + lage+ lexperience + lperformance,
  id = ~id,
  estimate = "mm")


df_lr_mm <- data.frame(lr_mm)
df_lr_mm<- df_lr_mm[,3:ncol(df_lr_mm)]

names(df_lr_mm) <- toTitleCase(tolower(names(df_lr_mm)))
df_lr_mm <- xtable(df_lr_mm) #latex table

print(df_lr_mm, type = "latex", include.rownames = FALSE, 
      caption.placement = "top")


lr_mm <- lr_mm[which(lr_mm$feature=="lorientation"),]
p2_lr <- ggplot(lr_mm,
                aes(x = estimate,
                    y = level)) +
  geom_point(position = position_dodge(width = 0.5),
             size = 2) +
  geom_errorbar(aes(xmin = lower, xmax = upper),
                position = position_dodge(width = 0.5),
                width = 0.2, size = 1) + 
  geom_vline(xintercept = 5) +
  xlim(4,6)+
  theme_bw() +
  theme(legend.position = "none") +
  ggtitle("MMs") +
  xlab(NULL) +
  labs(
    y = "",   
    x="Estimated MMs",
    color = "")

# 3x6
grid.arrange(p1_lr,p2_lr, ncol = 2)

# Standardized effects -----------------------------------------------------
rm(list = ls(all = TRUE))
load("final_data.RData")

final_data$subsleg_index_standardized <- scale(final_data$subsleg_index)

h1_amce <- cj(
  final_data,
  subsleg_index_standardized~ selectorate + candidates + polling + lsex + lorientation + lage + lexperience + lperformance,
  estimate = "amce",
  id = ~id
)

h1_amce <- subset(h1_amce, feature == "lorientation")
